#ifndef mesh_H
#define mesh_H

#include <list>
#include <vector>
#include <set>
#include <functional>
#include <math.h>
#include <algorithm>
#include <iostream>// must be deleted
#include <iomanip>// must be deleted

#include "excepts.h"
#include "geometry.h"
#include "datumhandler.h"

namespace ABAN
{
    class BoundaryContainer;

    #define CELL_CONTAINER_TYPE std::list
    #define CellIterator CELL_CONTAINER_TYPE<Cell>::iterator
    class Cell;
    class CellIteratorLessFunctional;


    //Cell Neighborhood Set
    typedef std::set<CellIterator,CellIteratorLessFunctional> CellNeighborhoodSet;
    #define CellNeighborhoodSetIterator CellNeighborhoodSet::iterator

    //a Cartesian Adaptive Cell
    struct Cell
    {
        //Data Fields
        SeqType no;                                     //cell number
        CellNeighborhoodSet* neighbors[facesPerCell];   //neighbors of cell in each face
        FloatType area[noOfAxis];                       //area of opposite faces are equal
        FloatType volume;                               //volume of cell
        Vector center;
        Vector lengths;
        FractionalCoordination fractCoord;
        SeqType grandParentNo;
        FloatType errMassConservation;
        Vector LSEgradP;
        Vector LSEgradV[noOfAxis];

        bool checked;

        //Methods
        Cell(); //default constructor
        Cell(const Cell& other);
        ~Cell();
        Cell& operator=(const Cell& other);

        void addThisToNeighbors(CellIterator thisCell);
        void deleteThisFromNeigbors(CellIterator thisCell);

        //functios for fractional coordination neighborhood.
        static bool overlap(FractionalCoordination &fc1,FractionalCoordination &fc2,SeqType dir);
        void correctInheritedNeighbors();
        bool isCoarsableWith(CellIterator aCell);
        //geometrical functions
        Vector getLSEGradient(DoubleDatumVector *fee);
        FloatType getNormalToFaceGradient(DoubleDatumVector *fee, SeqType face);
        FloatType getFaceValue(DoubleDatumVector *fee, SeqType face);

        Vector sharedFaceCenter(CellIterator aCell,SeqType dir);//assumed that these two cells are neghibors;
        FloatType sharedArea(CellIterator aCell,SeqType dir);//assumed that these two cells are neghibors;
        Vector getNode(SeqType n);
    };

    //Cell Container
    typedef CELL_CONTAINER_TYPE<Cell> RawCellContainer;
    class CellContainer : public RawCellContainer
    {
        //Note: According to number of data which each cell contains, a DatumContainer should be defined
        public:
        bool currCheckedVal,hasDirichletPBC;

        DatumContainer<CellConnectivity> *nodeItersVecContainer;
        CellContainer();
        ~CellContainer();

        void addMatchingFaceNeighbor(CellIterator cellIte1,CellIterator cellIte2, SeqType face);
        void checkForMatchingNeighbors();
        void freeConnectivityData();
        void appendGrandCell(SeqType cellNo,CellConnectivity &aGCell);
        //todo: adding a ghostcontainer manager class to manage all containers
        //it is needed specially for correcting u,v,w,p on cell break/coarse
        void breakCell(CellIterator iparent, SeqType dir);
        void coarseCell(CellIterator aCell, CellIterator bCell, SeqType dir);

        void updateGradients();

        DoubleDatumVector *vels[noOfAxis],*velsN[noOfAxis],*p,*pGradients[noOfAxis],;
        BoundaryContainer *boundary;
    };


    //CellIteratorLessFunctional : defines "<" operator for CellIterator
    class CellIteratorLessFunctional : public std::binary_function<CellIterator,CellIterator,bool>
    {
        public:
        result_type operator()(first_argument_type leftOperand,second_argument_type rightOperand);
    };

}//of namespace ABAN

#include "boundary.h"

//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////////
namespace ABAN
{
    /*******************************************************/
    /********************CellConnectivity*******************/
    /*******************************************************/
    void CellConnectivity::rearrangeNodes()
    {
        std::sort(begin(),end(),Vector::compare);
    }
    //-------------------------------------------------------
    bool CellConnectivity::isMatchingNeighborOnFace(CellConnectivity aGCell, SeqType face)
    {
        /*bool isNeighbor = true;
        SeqType oppFace = oppositeSide[face];

        //we used (nodesPerFace-1) because if only three nodes of two faces become equal,
        //those faces are neighbor and there is no need to check the last node
        for(int i = 0; (i<(nodesPerFace-1)) && isNeighbor;i++)
          isNeighbor = isNeighbor && (this->at(nodeNumbersInFaces[face   ][i]) ==
                                      aGCell.at(nodeNumbersInFaces[oppFace][i]));

        return isNeighbor;*/

        //even lesser condition is suffiicient: just equivalant of a pair of nodes on a particular position
        return (       at(nodeNumbersInFaces[face              ][0]) ==
                aGCell.at(nodeNumbersInFaces[oppositeSide[face]][0]));
    }
    /*******************************************************/
    /********************CellContainer**********************/
    /*******************************************************/
    CellContainer::CellContainer():currCheckedVal(true),hasDirichletPBC(false)
    {
        nodeItersVecContainer = new DatumContainer<CellConnectivity>();

        vels[dirXP]=new DoubleDatumVector();
        vels[dirYP]=new DoubleDatumVector();
        vels[dirZP]=new DoubleDatumVector();

        velsN[dirXP]=new DoubleDatumVector();
        velsN[dirYP]=new DoubleDatumVector();
        velsN[dirZP]=new DoubleDatumVector();

        p=new DoubleDatumVector();

        pGradients[dirXP]=new DoubleDatumVector();
        pGradients[dirYP]=new DoubleDatumVector();
        pGradients[dirZP]=new DoubleDatumVector();

        boundary=new BoundaryContainer();
    }
    //-------------------------------------------------------
    CellContainer::~CellContainer()
    {
          if (nodeItersVecContainer!=NULL) delete nodeItersVecContainer;
          delete vels[dirXP];delete vels[dirYP];delete vels[dirZP];
          delete velsN[dirXP];delete velsN[dirYP];delete velsN[dirZP];
          delete p;
          delete pGradients[dirXP];delete pGradients[dirYP];delete pGradients[dirZP];
          delete boundary;
    }
    //-------------------------------------------------------
    void CellContainer::checkForMatchingNeighbors()
    {
        SeqType face;
        bool matchedFound;

        CellIterator beginOfList = begin();
        CellIterator endOfList   = end();

        SeqType totalCells=size();
        SeqType processedCells=-1;

        cout<<setprecision(1);

        for (CellIterator toBeTested = beginOfList; toBeTested!=endOfList; toBeTested++)
        {
            for (face = dirXP; face<facesPerCell; face++)
            {
                for (CellIterator testWith = beginOfList; (!matchedFound) && (testWith!=endOfList); testWith++)
                    if (!nodeItersVecContainer->at(testWith->no).flag[oppositeSide[face]])
                        if (nodeItersVecContainer->at(testWith->no).isMatchingNeighborOnFace
                                    (nodeItersVecContainer->at(toBeTested->no),face))
                        {
                            matchedFound = true;
                            addMatchingFaceNeighbor(testWith,toBeTested,face);
                        }
                matchedFound = false;
            }// of face for
            if (++processedCells%10==0)
              cout<<"Processing neighborhoods... "<<(processedCells)/*/totalCells)*100<<"%"*/<<char(13)<<flush;
        }//of outer for
        //cout<<resetiosflags();
    }//of function
    //-------------------------------------------------------
    void CellContainer::addMatchingFaceNeighbor(CellIterator cellIte1,CellIterator cellIte2, SeqType face)
    {
        cellIte1->neighbors[face]->insert(cellIte2);
        nodeItersVecContainer->at(cellIte1->no).flag.set(face);

        cellIte2->neighbors[oppositeSide[face]]->insert(cellIte1);
        nodeItersVecContainer->at(cellIte2->no).flag.set(oppositeSide[face]);
    }
    //-------------------------------------------------------
    void CellContainer::freeConnectivityData()
    {
        delete nodeItersVecContainer;
        nodeItersVecContainer = NULL;
    }
    //-------------------------------------------------------
    void CellContainer::appendGrandCell(SeqType cellNo,CellConnectivity &aGCell)
    {
        Cell aCell;
        aCell.checked = currCheckedVal;
        aCell.no = cellNo;
        aCell.grandParentNo = cellNo;
        aCell.volume = 1.0;
        aCell.errMassConservation = 0.0;
        for (SeqType i=dirXP;i<noOfAxis;i++)
        {
            aCell.lengths[i] = (*aGCell[0])[i] - (*aGCell[7])[i];
            aCell.volume *= aCell.lengths[i];
            aCell.center[i] = ((*aGCell[0])[i] + (*aGCell[7])[i])/2.0;
        };
        aCell.area[dirXP] = aCell.lengths[dirYP] * aCell.lengths[dirZP];
        aCell.area[dirYP] = aCell.lengths[dirZP] * aCell.lengths[dirXP];
        aCell.area[dirZP] = aCell.lengths[dirXP] * aCell.lengths[dirYP];

        push_back(aCell);
        if (nodeItersVecContainer!=NULL) nodeItersVecContainer->push_back(aGCell);
    }
    //-------------------------------------------------------
    void CellContainer::updateGradients()
    {
        for(CellIterator cell=begin();cell!=end();cell++)
        {
            cell->LSEgradP = cell->getLSEGradient(p);
            for (SeqType dir=dirXP;dir<noOfAxis;dir++)
                cell->LSEgradV[dir] = cell->getLSEGradient(vels[dir]);
        }
    }
    //-------------------------------------------------------
    void CellContainer::breakCell(CellIterator iparent, SeqType dir)
    {
        iparent->volume /= 2.0;
        for (SeqType d=dirXP; d<noOfAxis; d++) iparent->area[d] /= 2.0;

        iparent->area[dir] *= 2.0;
        iparent->lengths[dir] /= 2.0;
        iparent->center[dir] -= iparent->lengths[dir] / 2.0;
        iparent->deleteThisFromNeigbors(iparent);

        Cell child(*iparent);
        child.no = (--end())->no + 1; //will be changed if node-numbering algorithms added to code
        child.center[dir] += child.lengths[dir];
        push_back(child);
        CellIterator ichild = --end();

        iparent->neighbors[dir              ]->clear();
        ichild ->neighbors[oppositeSide[dir]]->clear();
        iparent->neighbors[dir              ]->insert(ichild );
        ichild ->neighbors[oppositeSide[dir]]->insert(iparent);

        iparent->fractCoord[dir].split(true);
        ichild ->fractCoord[dir].split(false);

        iparent->correctInheritedNeighbors();
        ichild ->correctInheritedNeighbors();

        iparent->addThisToNeighbors(iparent);
        ichild ->addThisToNeighbors(ichild );
    }
    //-------------------------------------------------------
    void CellContainer::coarseCell(CellIterator aCell, CellIterator bCell, SeqType dir)
    {
        if (aCell->isCoarsableWith(bCell))
        {
            //Geometrical Correction
            aCell->volume *= 2.0;
            for (SeqType d=dirXP; d<noOfAxis; d++)  aCell->area[d] *= 2.0;
            aCell->area[dir] /= 2.0;
            aCell->lengths[dir] *= 2.0;
            aCell->center[dir] = (aCell->center[dir]+bCell->center[dir]) / 2.0;

            //Neighborhood Correction
            bCell->deleteThisFromNeigbors(bCell);
            aCell->neighbors[dir              ]->erase(bCell);
            bCell->neighbors[oppositeSide[dir]]->erase(aCell);
            bCell->addThisToNeighbors(aCell);

            for (SeqType d=dirXP;d<noOfDirections;d++)
              aCell->neighbors[d]->insert(bCell->neighbors[d]->begin(),bCell->neighbors[d]->end());

            //Renumber

            aCell->fractCoord[dir].collapse();

            erase(bCell);
        }
    }
    /*******************************************************/
    /***********************Cell****************************/
    /*******************************************************/
    Cell::Cell()
    {
        for (SeqType i=dirXP;i<facesPerCell;i++) neighbors[i] = new CellNeighborhoodSet();
    }
    //-------------------------------------------------------
    Cell::~Cell()
    {
        for (SeqType i=dirXP;i<facesPerCell;i++) delete neighbors[i];
    }
    //-------------------------------------------------------
    Cell::Cell(const Cell& other)
    {
        for (SeqType i=dirXP;i<facesPerCell;i++) neighbors[i] = new CellNeighborhoodSet();
        *this = other;
    };
    //-------------------------------------------------------
    Cell& Cell::operator=(const Cell& other)
    {
        no=other.no;
        volume=other.volume;
        center=other.center;
        lengths=other.lengths;
        grandParentNo=other.grandParentNo;
        checked=other.checked;
        LSEgradP=other.LSEgradP;

        for (SeqType i=dirXP;i<facesPerCell;i++) *neighbors[i]=*other.neighbors[i];
        for (SeqType i=dirXP;i<noOfAxis;i++)
        {
            LSEgradV[i]=other.LSEgradV[i];
            area[i]=other.area[i];
            fractCoord[i]=other.fractCoord[i];
        }

        return *this;
    }
    //-------------------------------------------------------
    void Cell::addThisToNeighbors(CellIterator thisCell)
    {
        for (SeqType face = dirXP; face<facesPerCell; face++)
          for (CellNeighborhoodSetIterator neighbor = neighbors[face]->begin();
               neighbor!=neighbors[face]->end();
               neighbor++)

              (*neighbor)->neighbors[oppositeSide[face]]->insert(thisCell);
    }
    //-------------------------------------------------------
    void Cell::deleteThisFromNeigbors(CellIterator thisCell)
    {
        for (SeqType face = dirXP; face<facesPerCell; face++)
            for (CellNeighborhoodSetIterator neighbor = neighbors[face]->begin();
                 neighbor!=neighbors[face]->end();
                 neighbor++)

            (*neighbor)->neighbors[oppositeSide[face]]->erase(thisCell);
    }
    //-------------------------------------------------------
    bool Cell::overlap(FractionalCoordination &fc1,FractionalCoordination &fc2,SeqType dir)
    {
        bool result = true;
        for (int i = dirXP;i<noOfAxis;i++)
            if (i!=dir) result = (result && (fc1[i].relation(fc2[i])!=fcrNone));
        return result;
    }
    //-------------------------------------------------------
    void Cell::correctInheritedNeighbors()
    {
      for (SeqType face = dirXP; face<facesPerCell; face++)
          for (CellNeighborhoodSetIterator neighb = neighbors[face]->begin();
               neighb != neighbors[face]->end();
               neighb++)
              if (!overlap(this->fractCoord,
                           (*neighb)->fractCoord,
                           face<noOfAxis?face:oppositeSide[face]))

                  neighbors[face]->erase(neighb);
    }
    //-------------------------------------------------------
    bool Cell::isCoarsableWith(CellIterator aCell)
    {
        bool result=(grandParentNo==aCell->grandParentNo);

        if (result)
            for (SeqType dir = 0; dir<noOfAxis; dir++)
                result = result && (fractCoord[dir].denominator == aCell->fractCoord[dir].denominator);

        return result;
    }
    //-------------------------------------------------------
    Vector Cell::getLSEGradient(DoubleDatumVector *fee)
    {
        FloatType x2bar=0.0,xybar=0.0,xzbar=0.0,xbar=0.0
                           ,y2bar=0.0,yzbar=0.0,ybar=0.0
                                     ,z2bar=0.0,zbar=0.0;
        FloatType xfeebar=0.0,yfeebar=0.0,zfeebar=0.0,feebar=0.0;
        FloatType xbar2=0.0,ybar2=0.0,zbar2=0.0,
                  xzbar2=0.0, xybar2=0.0, yzbar2=0.0,
                  denominator=0.0;
        bool hasNeighbOnDir[noOfAxis] = {false, false, false};
        SeqType neighbNum=1;//1 because the current cell itself included in average computation

        Vector result;

        for (SeqType face = dirXP; face<facesPerCell; face++)
            for (CellNeighborhoodSetIterator neighb = neighbors[face]->begin();
                 neighb != neighbors[face]->end();
                 neighb++)
            {
                hasNeighbOnDir[(face<noOfAxis?face:oppositeSide[face])] = true;
                neighbNum++;
                SeqType nno=(*neighb)->no;

                xbar+=(*neighb)->center[dirXP];
                ybar+=(*neighb)->center[dirYP];
                zbar+=(*neighb)->center[dirZP];

                x2bar+=sqr((*neighb)->center[dirXP]);
                y2bar+=sqr((*neighb)->center[dirYP]);
                z2bar+=sqr((*neighb)->center[dirZP]);

                xybar+=(*neighb)->center[dirXP]*(*neighb)->center[dirYP];
                xzbar+=(*neighb)->center[dirXP]*(*neighb)->center[dirZP];
                yzbar+=(*neighb)->center[dirYP]*(*neighb)->center[dirZP];

                feebar+=(*fee)[nno];
                xfeebar+=(*neighb)->center[dirXP]*(*fee)[nno];
                yfeebar+=(*neighb)->center[dirYP]*(*fee)[nno];
                zfeebar+=(*neighb)->center[dirZP]*(*fee)[nno];
            }

        xbar=(center[dirXP]+xbar)/neighbNum;
        ybar=(center[dirYP]+ybar)/neighbNum;
        zbar=(center[dirZP]+zbar)/neighbNum;

        x2bar=(x2bar+sqr(center[dirXP]))/neighbNum;
        y2bar=(y2bar+sqr(center[dirYP]))/neighbNum;
        z2bar=(z2bar+sqr(center[dirZP]))/neighbNum;

        xybar=(xybar+center[dirXP]*center[dirYP])/neighbNum;
        xzbar=(xzbar+center[dirXP]*center[dirZP])/neighbNum;
        yzbar=(yzbar+center[dirYP]*center[dirZP])/neighbNum;

        feebar=(feebar+(*fee)[no])/neighbNum;
        xfeebar=(xfeebar+center[dirXP]*(*fee)[no])/neighbNum;
        yfeebar=(yfeebar+center[dirYP]*(*fee)[no])/neighbNum;
        zfeebar=(zfeebar+center[dirZP]*(*fee)[no])/neighbNum;

        if (!hasNeighbOnDir[dirXP]){x2bar=1.0;xfeebar=xbar=xzbar=xybar=0.0;}
        if (!hasNeighbOnDir[dirYP]){y2bar=1.0;yfeebar=ybar=xybar=yzbar=0.0;}
        if (!hasNeighbOnDir[dirZP]){z2bar=1.0;zfeebar=zbar=xzbar=yzbar=0.0;}

        xbar2=sqr(xbar);
        ybar2=sqr(ybar);
        zbar2=sqr(zbar);

        xzbar2=sqr(xzbar);
        xybar2=sqr(xybar);
        yzbar2=sqr(yzbar);

        denominator=(x2bar*y2bar*z2bar - x2bar*y2bar*zbar2 - x2bar*yzbar2 + 2.0*x2bar*yzbar*ybar*zbar - x2bar*ybar2*z2bar - xybar2*z2bar + xybar2*zbar2 + 2.0*xybar*yzbar*xzbar - 2.0*xybar*yzbar*xbar*zbar - 2.0*xybar*ybar*xzbar*zbar + 2.0*xybar*ybar*xbar*z2bar - y2bar*xzbar2 + 2.0*xzbar*y2bar*xbar*zbar + xzbar2*ybar2 - 2.0*xzbar*ybar*xbar*yzbar - y2bar*xbar2*z2bar + xbar2*yzbar2);

        result[dirXP] = - ( - xfeebar*y2bar*z2bar + xfeebar*y2bar*zbar2 + xfeebar*yzbar2 - 2.0*xfeebar*yzbar*ybar*zbar + xfeebar*ybar2*z2bar + yfeebar*xybar*z2bar - yfeebar*xybar*zbar2 - yfeebar*yzbar*xzbar + yfeebar*yzbar*xbar*zbar + yfeebar*ybar*xzbar*zbar - yfeebar*ybar*xbar*z2bar - zfeebar*xybar*yzbar + zfeebar*xybar*ybar*zbar + zfeebar*y2bar*xzbar - zfeebar*y2bar*xbar*zbar - zfeebar*xzbar*ybar2 + zfeebar*ybar*xbar*yzbar + feebar*xybar*yzbar*zbar - feebar*xybar*ybar*z2bar - feebar*y2bar*xzbar*zbar + feebar*y2bar*xbar*z2bar + feebar*yzbar*xzbar*ybar - feebar*xbar*yzbar2) / denominator;
        result[dirYP] = - (   xfeebar*xybar*z2bar - xfeebar*xybar*zbar2 - xfeebar*yzbar*xzbar + xfeebar*yzbar*xbar*zbar + xfeebar*ybar*xzbar*zbar - xfeebar*ybar*xbar*z2bar - yfeebar*x2bar*z2bar + yfeebar*x2bar*zbar2 + yfeebar*xzbar2 - 2.0*yfeebar*xzbar*xbar*zbar + yfeebar*xbar2*z2bar + zfeebar*x2bar*yzbar - zfeebar*x2bar*ybar*zbar - zfeebar*xzbar*xybar + zfeebar*xzbar*xbar*ybar + zfeebar*xbar*xybar*zbar - zfeebar*xbar2*yzbar - feebar*x2bar*yzbar*zbar + feebar*x2bar*ybar*z2bar + feebar*xzbar*xybar*zbar - feebar*xzbar2*ybar - feebar*xbar*xybar*z2bar + feebar*xbar*yzbar*xzbar) / denominator;
        result[dirZP] = - ( - xfeebar*xybar*yzbar + xfeebar*xybar*ybar*zbar + xfeebar*y2bar*xzbar - xfeebar*y2bar*xbar*zbar - xfeebar*xzbar*ybar2 + xfeebar*ybar*xbar*yzbar + yfeebar*x2bar*yzbar - yfeebar*x2bar*ybar*zbar - yfeebar*xzbar*xybar + yfeebar*xzbar*xbar*ybar + yfeebar*xbar*xybar*zbar - yfeebar*xbar2*yzbar - zfeebar*x2bar*y2bar + zfeebar*x2bar*ybar2 + zfeebar*xybar2 - 2.0*zfeebar*xybar*xbar*ybar + zfeebar*xbar2*y2bar + feebar*x2bar*y2bar*zbar - feebar*x2bar*ybar*yzbar - feebar*xybar2*zbar + feebar*xybar*xzbar*ybar + feebar*xbar*xybar*yzbar - feebar*xbar*y2bar*xzbar) / denominator;

        return result;
    }
    //-------------------------------------------------------
    FloatType Cell::getNormalToFaceGradient(DoubleDatumVector *fee, SeqType face)
    {
        SeqType pCellNo=no;
        Vector rP=center;
        FloatType result=0.0;
        Vector gradientFeeP = getLSEGradient(fee);
        SeqType dir = (face<noOfAxis?face:oppositeSide[face]);
        for (CellNeighborhoodSetIterator neighbor=neighbors[face]->begin();
             neighbor!=neighbors[face]->end();
             neighbor++)
        {
            SeqType nCellNo = (*neighbor)->no;
            Vector rN = (*neighbor)->center;
            FloatType deltaRDir = (rN[dir]-rP[dir]);
            Vector rF = sharedFaceCenter((*neighbor),dir);
            Vector gradientFeeN=(*neighbor)->getLSEGradient(fee);
            FloatType areaF = sharedArea((*neighbor),dir);
            FloatType addingPart = (*fee)[nCellNo] - (*fee)[pCellNo];
            for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                if (otherDir!=dir)
                    addingPart+=gradientFeeN[otherDir]*(rF[otherDir]-rN[otherDir])-
                                gradientFeeP[otherDir]*(rF[otherDir]-rP[otherDir]);
            result+=addingPart/fabs(deltaRDir)*areaF;
        }
        return (result/area[dir]);
    }
    //-------------------------------------------------------
    FloatType Cell::getFaceValue(DoubleDatumVector *fee, SeqType face)
    {
        SeqType pCellNo=no;
        Vector rP=center;
        FloatType result=0.0;
        Vector gradientFeeP = getLSEGradient(fee);
        SeqType dir = (face<noOfAxis?face:oppositeSide[face]);
        for (CellNeighborhoodSetIterator neighbor=neighbors[face]->begin();
             neighbor!=neighbors[face]->end();
             neighbor++)
        {
            SeqType nCellNo = (*neighbor)->no;
            Vector rN = (*neighbor)->center;
            FloatType deltaRDir = (rN[dir]-rP[dir]);
            Vector rF = sharedFaceCenter((*neighbor),dir);
            Vector gradientFeeN = (*neighbor)->getLSEGradient(fee);
            FloatType areaF = sharedArea((*neighbor),dir);
            FloatType feeF=(rF[dir]-rP[dir])*(*fee)[nCellNo]+(rN[dir]-rF[dir])*(*fee)[pCellNo];
            for(SeqType otherDir=dirXP;(otherDir<noOfAxis);otherDir++)
                if (otherDir!=dir)
                    feeF+=(rF[dir]-rP[dir])*(rF[otherDir]-rN[otherDir])*gradientFeeN[otherDir]+
                          (rN[dir]-rF[dir])*(rF[otherDir]-rP[otherDir])*gradientFeeP[otherDir];

            result+=feeF/deltaRDir*areaF;
        }
        return (result/area[dir]);
    }
    //-------------------------------------------------------
    Vector Cell::sharedFaceCenter(CellIterator aCell,SeqType dir)
    {
        Vector result;

        for(SeqType adir=dirXP;adir<noOfAxis;adir++)
            if (adir == dir)
            {
                FloatType coeff = (center[adir]>aCell->center[adir]?-1.0:1.0);
                result[adir]=center[adir]+coeff*lengths[adir]/2.0;
            }
            else
                if (fractCoord[adir].denominator<aCell->fractCoord[adir].denominator)
                    result[adir]=center[adir];
                else
                    result[adir]=aCell->center[adir];
        return result;
    }
    //-------------------------------------------------------
    FloatType Cell::sharedArea(CellIterator aCell,SeqType dir)
    {
      SeqType s1,b2=0;
      if (dir!=dirXP) s1 = fractCoord[dirXP].denominator; else s1 = fractCoord[dirYP].denominator;
      for(SeqType adir=dirXP;adir<noOfAxis;adir++)
          if (adir != dir)
          {
              if (fractCoord[adir].denominator<s1) s1 = fractCoord[adir].denominator;
              if (aCell->fractCoord[adir].denominator>b2) b2 = fractCoord[adir].denominator;
          }
      return static_cast<FloatType>(s1)/static_cast<FloatType>(b2)*area[dir];
    }
    //-------------------------------------------------------
    Vector Cell::getNode(SeqType n)
    {
        Vector result;
        result=center;
        for (SeqType i=dirXP; i<noOfAxis; i++)
            result[i]+=(facesOfNodes[n][i]<noOfAxis?lengths[i]/2.0:-lengths[i]/2.0);
        return result;
    }
    /*******************************************************/
    /***************CellIteratorLessFunctional**************/
    /*******************************************************/
    CellIteratorLessFunctional::result_type CellIteratorLessFunctional::operator()
      (CellIteratorLessFunctional::first_argument_type leftOperand,
       CellIteratorLessFunctional::second_argument_type rightOperand)
    {
        result_type result;
        FloatType R;

        R = (leftOperand->center[dirXP] - rightOperand->center[dirXP]);result = R<0;
        if (fabs(R) < zeroRadious)
        {
            R = (leftOperand->center[dirYP] - rightOperand->center[dirYP]);result = R<0;
            if (fabs(R) < zeroRadious)
            {
                R = (leftOperand->center[dirZP] - rightOperand->center[dirZP]);
                result = (R<0) && !(fabs(R)<zeroRadious);
            }
        }

        return result;
    }
    /*******************************************************/

}


#endif
